BiocManager::install(c("DESeq2", "SummarizedExperiment", "S4Vectors", "BiocGenerics", "IRanges", "GenomicRanges"), force = TRUE, update = TRUE)
library(DESeq2)
library(ggplot2)

setwd("D:/1大学生活/大三下/zxy/R语言/class2")
# 读取临床信息数据
clin_inf <- read.table("clin_inf.csv", header = TRUE)

# 读取基因表达计数数据
count_data <- read.table("count.csv", header = TRUE, row.names = 1, check.names = FALSE)

# 将数据转换为整数
count_data <- round(count_data)
count_data <- as.data.frame(apply(count_data, 2, as.integer))
rownames(count_data) <- rownames(read.csv("count.csv", header = TRUE, row.names = 1))

# 数据预处理和整理
# 1. 提取样本名并匹配临床信息
sample_names <- colnames(count_data)
# 处理样本名，移除末尾的字母（如"D"、"B"、"A"）
clean_sample_names <- str_replace(sample_names, "[A-Z]$", "")

# 2. 创建实验设计矩阵(metadata)
metadata <- data.frame(sample = sample_names, stringsAsFactors = FALSE)
metadata$clean_name <- clean_sample_names

# 3. 与临床信息合并
metadata <- merge(metadata, clin_inf, by.x = "clean_name", by.y = "样本名称", all.x = TRUE)

# 4. 根据年龄分组（65岁以前为中年组，65岁以后为老年组）
metadata$age_group <- ifelse(metadata$年龄 < 65, "中年组", "老年组")
metadata$age_group <- factor(metadata$age_group)

# 5. 确保所有样本都有对应的分组信息
metadata <- metadata[!is.na(metadata$age_group), ]

# 6. 仅保留有完整信息的样本
valid_samples <- metadata$sample
count_data_filtered <- count_data[, valid_samples]

# 7. 确保metadata行顺序与count_data列顺序一致
metadata <- metadata[match(colnames(count_data_filtered), metadata$sample), ]
rownames(metadata) <- metadata$sample

# 8. 构建DESeq2对象进行差异分析
dds <- DESeqDataSetFromMatrix(
  countData = count_data_filtered,
  colData = metadata,
  design = ~ age_group
)


# 10. 设置因子水平，将"中年组"设为参照组
dds$age_group <- relevel(factor(dds$age_group), ref = "中年组")

# 11. 运行DESeq2分析
dds <- DESeq(dds)

# 12. 提取差异分析结果
res <- results(dds, contrast = c("age_group", "老年组", "中年组"))

res_df <- as.data.frame(res)

# 13. 添加基因ID作为列
res_df$gene_id <- rownames(res_df)

# 16. 数据可视化
# 添加分类标签
res_df$diffexpressed <- "NO"
res_df$diffexpressed[res_df$padj < 0.05 & res_df$log2FoldChange > 1 & !is.na(res_df$padj)] <- "UP"
res_df$diffexpressed[res_df$padj < 0.05 & res_df$log2FoldChange < -1 & !is.na(res_df$padj)] <- "DOWN"

# 过滤无效数据（NA和Inf）
res_df_filtered <- res_df %>%
  filter(
    !is.na(padj),
    is.finite(-log10(padj))  # 移除padj=0导致的Inf值
  )

# 绘制火山图（y轴从0开始）
volcano_plot <- ggplot(res_df_filtered, aes(x = log2FoldChange, y = -log10(padj))) +
  geom_point(aes(color = diffexpressed), alpha = 0.6, size = 1.5) +
  scale_color_manual(
    values = c("UP" = "red", "DOWN" = "blue", "NO" = "grey"),
    labels = c("Downregulated", "Not significant", "Upregulated")
  ) +
  geom_vline(xintercept = c(-1, 1), linetype = "dashed", color = "black", linewidth = 0.5) +
  geom_hline(yintercept = -log10(0.05), linetype = "dashed", color = "black", linewidth = 0.5) +
  scale_y_continuous(limits = c(0.1, 5)) +  # 设置明确上限，比如6
  labs(
    x = expression(log[2]("Fold Change")),
    y = expression(-log[10]("Adjusted p-value")),
    title = "Volcano Plot: 老年组 vs 中年组",
    color = "Differential Expression"
  ) +
  theme_classic() +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold", size = 14),
    axis.title = element_text(size = 12),
    axis.text = element_text(size = 10),
    legend.position = "right",
    legend.title = element_text(size = 10),
    legend.text = element_text(size = 8)
  )

ggsave("volcano_plot.png", plot = volcano_plot, width = 8, height = 6, dpi = 300)

